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Abstract 

Given a set S = {vi, . . . ,Vn} C and a point p £ R™, testing if p £ conv{S), the convex hull 
of S, is a fundamental problem in computational geometry and linear programming. First, we prove a 
Euclidean distance duality: Either for each p' € conv{S) \ {p}, 3 Vj £ S such that d{p' ,Vj) > d{p,Vj); or 
3 p' £ conv{S) such that d{p',Vi) < d{p,Vi),yi. The former case holds if and only if p € conv{S) and 
the latter if and only if p conv{S). Utilizing this duality, we describe a simple fully polynomial time 
approximation scheme, called the Triangle Algorithm: Given e > 0, in 0{mne~^) arithmetic operations 
it computes p' G conv(S) where, either d{p',p) < ed{p,Vj) for some j; or d{p',Vi) < d{p,Vi),yi. When 
p ^ conv{S), we show d{p,p') estimates the distance from p to conv{S) to within a factor of two. We also 
show how to solve general LP via the Triangle Algorithm and give a corresponding complexity analysis. 
In particular, we prove a sensitivity theorem that converts LP feasibility with bounded domain into a 
convex hull problem, then gives the necessary accuracy for computing an e-approximate solution. We 
also contrast the theoretical performance of the Triangle Algorithm with the sparse greedy approximation 
(equivalent to Frank- Wolfe and Gilbert algorithms) for the minimization of a convex quadratic over a 
simplex, a problem arising in machine learning, approximation theory, and statistics. 

Keywords: Convex Hull, Linear Programming, Duality, Approximation Algorithms, Sparse Greedy 
Approximation, Frank- Wolfe Algorithm, Support Vector Machines. 

1 Introduction 

Given a set S = {vi, ...,«„} C M™ and a distinguished point p G M™, we consider the problem of testing if 
p lies in conv{S), the convex hull of 5*. Throughout the article we shall refer to this problem as the convex 
hull decision problem, or simply as problem (P). This is a basic problem in computational geometry and a 
very special case of the convex hull problem, a problem that according to Goodman and O'Rourke |17| . is a 
"catch-all phrase for computing various descriptions of a polytope that is either specified as the convex hull 
of a finite point set or the intersection of a finite number of halfspaces." The descriptions include those of 
vertices, facets, and adjacencies. 

Problem (P) is not only a fundamental problem in computational geometry, but in linear programming 
(LP). This can be argued on different grounds. On the one hand problem (P) is a very special case of LP. 
On the other hand, it is well known that through LP duality theory the general LP problem may be cast as 
a single LP feasibility problem, see e.g. Chvatal |7|. The LP feasibility problem can then be converted into 
problem (P) via several different approaches. To argue the significance of (P) in linear programming, it can 
be justified that the two most famous polynomial-time LP algorithms, the ellipsoid algorithm of Khachiyan 
P7] and the projective algorithm of Karmarkar J^S] , are in fact explicitly or implicitly designed to solve a 
case of problem (P) where p — see Qjj. Furthermore, using an approach suggested by Chvatal, in [H] 
it is shown that there is a direct connection between a general LP feasibility and this homogeneous case of 
problem (P), over integer or real input data. For integer inputs all known polynomial-time LP algorithms 
- when applied to solve problem (P) - would exhibit a theoretical complexity that is polynomial in m, n, 
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and the size of encoding of the input data, often denoted by L, see e.g. [57]. The number L can generally 
be taken to be dependent on the logarithm of mn and of the logarithm of the absolute value of the largest 
entry in the input data. These are sometimes known as weakly polynomial-time algorithms. No strongly 
polynomial-time algorithm is known for LP, i.e. an algorithm that would in particular solve problem (P) in 
time complexity polynomial in m and n. 

Problem (P) finds applications in computational geometry itself, in particular among the variations of 
the convex hull problem. One such a variation is the irredundancy problem, the problem of computing 
all the vertices of conv{S), see [T7]. Clearly, any algorithm for LP can be used to solve the irredundancy 
problem by solving a sequence of 0(n) convex hull decision problems. Clarkson [3], has given a more efficient 
algorithm than the straightforward approach of solving n linear programs. The complexity of his algorithm 
depends on the number of vertices of the convex hull. Furthermore, by using an approach originally due to 
Matousek [311 for solving closely related linear programming problems, Chan ^ has given a faster algorithm 
for the irredundancy problem. What is generally considered to be the convex hull problem is a problem 
more complicated than (P) and the irredundancy problem, requiring the description of conv{S) in terms of 
its vertices, facets and adjacencies, see Chazelle [6 who offers an optimal algorithm. 

Problem (P) can also be formulated as the minimization of a convex quadratic function over a simplex. 
This particular convex program has found applications in statistics, approximation theory, and machine 
learning, see e.g Clarkson [10] and Zhang 39_ who consider the analysis of a greedy algorithm for the more 
general problem of minimizing certain convex functions over a simplex (equivalently, maximizing concave 
functions over a simplex). The oldest version of such greedy algorithms is Frank- Wolfe algorithm, [H] . 
Special cases of the problem include support vector machines (SVM), approximating functions as convex 
combinations of other functions, see e.g. Clarkson [TU]. The problem of computing the closest point of the 
convex hull of a set of points to the origin, known as polytope distance is the case of problem (P) where p is the 
origin. In some applications the polytope distance refers to the distance between two convex hulls. Gilbert's 
algorithm [16j for the polytope distance problem is one of the earliest known algorithms for the problem. 
Gartner and Jaggi [15] show Gilbert's algorithm coincides with Frank- Wolfe algorithm when applied to the 
minimization of a convex quadratic function over a simplex. In this case the algorithm is known as sparse 
greedy approximation. For many other results regarding the applications of the minimization of a quadratic 
function over a simplex, see the bibliographies in [39], [10] and [l^. Clarkson [10] gives a through analysis 
of Frank- Wolfe and its variations. 

Despite all the previous approaches and achievements in solving problem (P), in various formulations and 
via various algorithms, investigation of the problem will most likely continue in the future. The present article 
focuses on solving problem (P) directly. In the process it reveal new and interesting geometric properties 
of this fundamental convex hull problem, including a characterization theorem, leading to a new separating 
hyperplane theorem for problem (P), a distance duality theorem. It then utilizes the distance duality to 
describe a very simple algorithm called. Triangle Algorithm. Next, it utilizes the Triangle Algorithm, together 
with a sensitivity theorem, to give a new algorithm for the LP feasibility problem when there are no recession 
directions. It then analyzes the complexity of solving the general LP feasibility via the Triangle Algorithm. 

Before formally describing and proving the results in Section[2]- Section[6l in the remaining of this section 
we give an overview of the results. In 11.11 we describe the Triangle Algorithm and its properties. In 11.21 
we describe the geometry of the Triangle Algorithm and two geometric problems that can be associate as 
problems dual to problem (P). lu ll. 31 we review the literature on the properties and complexity of a greedy 
algorithm for the minimization of convex quadratic over a simplex and contrast these with the theoretical 
performance of the Triangle Algorithm. In 11.41 we describe the application of the Triangle Algorithm for 
solving an LP feasibility problem having no recession direction. In 11.51 we describe connections between 
problem (P) and general LP feasibility and LP optimization, as well as reviewing some of the vast literature 
on LP. lu ll. 61 we give the outline of the formal results in the remaining sections. 

1.1 The Triangle Algorithm and Its Properties 

Denoting the Euclidean distance between u, w G M™ by d{u, v) — '\/X]"=i('"« ~ the Triangle Algorithm 
takes as input a set of point S — {vi, . . . , w„} and a point p in R™. It consists of iterating two steps: 
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Triangle Algorithm {S = {vi, . . . , u„}, p) 

• Step 1. Given p' G conv{S) \ {p}, check if tliere exists Vj G S sucli tliat d{p' , Vj) > d{p, Vj). 
If no such Vj exists, stop, p ^ conv{S). 

• Step 2. Otherwise, on the line segment joining p' to Vj compute the point nearest to p. 
Denote this by p" . Replace p' with p", go to Step 1. 



Clearly, if the algorithm does not terminate in Step 1 it will loop indefinitely. We prove: Given e G 
(0,1), the Triangle Algorithm in at most 48mne^^ = 0{mne^^) arithmetic operations computes a point 
p' G conv{S) such that either 

d{p',p) < ed{p, Vj), for some j, (1) 

or 

d{p',Vi) < d{p,Vi), Vi = l,...,n. (2) 

Remark 1. Note that the constant in the worst-case complexity is an absolute constant and is independent 
of the input data. This is due to the fact that the Triangle Algorithm seeks to reduce the relative error, 
d{p' ,p)/d{p, Vj), which is a more sensible measure for the problem than seeking to reduce the absolute error, 
d{p' ,p). Clearly, approximation to a prescribed absolute error can also be achieved. 

Remark 2. By squaring the distances we have 

d{p', V,) > dip, V,) ^ dip', 0)2 - dip, 0)2 > 2vfip' - p). (3) 

Thus Step 1 does not require taking square-roots. Also, the computation of p" requires no square-root 
operations. 

We refer to a point p' satisfying ([T|) as iterate and Vj as pivot point. We refer to a point p' satisfying ([2]) 
as witness. This condition holds if and only ii p ^ conviS). This is because in this case we can prove the 
Voronoi cell of p' with respect to the two point set {p,p'} contains conviS) (see Figure[T]). Equivalently, the 
orthogonal bisector of the line segment pp' separates p from conviS). 

The set Wp of all such witnesses is the intersection of conviS) and the open balls, Bi ^ {x E M™ : 
dix, Vi) < Vi}, i ^ I, . . . ,n. Wp is a convex open set in the relative interior of conviS) (see Figure[3]). 

The justification in the name of the algorithm lies in the fact that in each iteration the algorithm searches 
for a triangle App'vj where Vj € S, p' E conviS) \ {p}, such that dip' , Vj) > dip, Vj). Given that such triangle 
exists, it uses Vj as a pivot point to "pull" the current iterate p' closer to p to get a new iterate p" G conviS). 

The correctness of the Triangle Algorithm lies in a new duality (theorem of the alternative) we call 
distance duality: 

Distance Duality 

Precisely one of the two conditions is satisfied: 

(i) : For each p' G conviS) \ {p}, there exists Vj G S such that dip' , Vj) > dip, Vj); 

(ii) : There exists p' G conviS) such that dip' , Vi) < dip, Vi), for alH = 1, . . . , n. 



The first condition is valid if and only if p G conviS), and the second condition if and only ifp ^ conviS). 
From the description of the Triangle Algorithm we see that given a point p' G conviS) that is not a witness, 
having dip,p') as the current gap, the Triangle Algorithm moves to a new point p" G conviS) where the 
new gap dip,p") is reduced. We will prove that when p G conv{S), the number of iterations K^, needed to 
get an approximate solution p' satisfying ([T]) is bounded above by 48e~^ = Oie~^). In the worst-case each 
iteration of Step 1 requires 0(mn) arithmetic operations. However, it may also take only 0(m) operations. 
The number of arithmetic operations in each iteration of Step 2 is only 0(m). Thus the Triangle Algorithm 
is a fully polynomial-time approximation scheme whose complexity for computing an e-approximate solution 
is 0(mne~2) arithmetic operation. In particular, for fixed e the complexity of the algorithm is only 0(mn). 
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The worst-case iteration complexity estimate is under the assumption of worst-case performance in each 
iteration of the algorithm. Hence, in practice a more efhcient complexity is to be expected. The algorithm 
is geometric in nature, simple in description, and very easy to implement. 

When p ^ conv{S), the Triangle Algorithm does not attempt to compute the closest point to p, say 
p„ G conv{S), rather a separating hyperplane. However, by virtue of the fact that it finds a special separating 
hyperplane, i.e. a hyperplane orthogonally bisecting the line pp' , it in the process computes an approximation 
to d{p,p^:) to within a factor of two. More precisely, any witness p' satisfies the inequality 

.bd{p,p')<d{p,p^)<d{p,p'). (4) 

Remark 3. Not only this approximation is useful for problem (P), but the case of computing the distance 
between two convex hulls, the polytope distance problem. As is well known the Minkowski difference of two 
convex hulls is a polytope whose shortest vector has norm equal to the distance between the two polytopes, 
see e.g. Clarkson [10] and Gartner and Jaggi [15j . 

1.2 The Geometry of the Triangle Algorithm 

To arrive at the distance duality mentioned above, we first prove a characterization theorem that leads to 
this new theorem of the alternative for problem (P). We remark here that the distance duality theorem is 
distinct from the classical Farkas lemma, or Gordan theorem. To arrive at this theorem we first prove: 

p e conv{S) if and only if given any point p' G conv{S) \ {p}, there exists vj such that d{p', vj) > d{p, vj). 

Next, we show the strict inequality can be replaced with d{p' ,Vj) > d{p,Vj). Thus the contrapositive 
theorem is: 

p ^ conv{S) if and only if there exists p' £ conv{S) such that d(p' , Vi) < d{p, Vi), for all i — 1, . . . ,n. 

These two results together imply the distance duality. A corollary of our characterization theorem reveals 
a geometric property of a set of balls. To describe this property, denote an open ball B, its closure B, and 
its boundary by dB, i.e. 

B = {.T e M" : d{x, v) < r}, B ^ {x € M™ : d{x, v) < r}, dB ^ {x e M" : d{x, v) = r}. 
Consider a set of open balls Bi — {x £ R™ : d(x, Vi) < 7'i}, i = 1, . . . ,n, and let S = {vi, . . . , ?;„}. 



Intersecting Balls Property: 




If p e nf^j^dB,, then 




p e conv{S) <=^ (n-LiSj) n conv{S) = < 


=^ {r\'}^iB,)r\conv{S) = {p}. 



In words, suppose a set of open balls have a common boundary point p. Then p lies in the convex hull 
of their centers, if and only if the intersection of the open balls is empty, if and only ii p is the only point in 
the intersection of the closure of the balls. A depiction of this property for a triangle is given in Figure [21 
This property suggests we can define a geometric dual for problem (P): 

Problem (Q) (Intersecting Balls Problem): 

Suppose there exists p e C\^=idBi. Determine if (n^^^Bi) fl conv{S) is nonempty. 



In fact the intersecting balls problem can be stated in more generality: 

Problem (Q') (General Intersecting Balls Problem): 

Suppose there exists p G ri'l^^dBi. Determine if (n"^^^^) i*^ nonempty. 



The Triangle Algorithm results in a fully polynomial-time approximation scheme for solving problems 
(Q) and (Q') with the same time complexity as that for problem (P). 
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When a point p lies in conv{S) fl (n"^;^9i?i), the union of the balls, U"^ji3i is referred as the forbidden 
zone of the convex hull of the centers, see [11] or [12] for the definition and some of its properties. The notion 
of forbidden zone of a convex set is significant and intrinsic in the characterization of the so-called mollified 
zone diagrams^ a variation of zone diagram of a finite set of points in the Euclidean plane, see [12) . The 
notion of zone diagram, introduced by Asano et al is itself a very rich and interesting variation of the 
classical Voronoi diagram, see e.g. [5], [Mj. Forbidden zones help give a characterization of mollified zone 
diagrams, in particular a zone diagram, |12j . For some geometric properties of forbidden zones of polygons 
and polytopes, see [1]. 

1.3 Comparison of Triangle Algorithm and Sparse Greedy Approximation 

Formally, the distance between p and conv{S) is defined as 

A = min{d{p',p) : p' E conv{S)} — d{p^,,p)}. (5) 

We have, p ^ conv{S), if and only if A > 0. While problem (P) does not require the computation 
of A when it is positive, in some applications this distance is required. However, as stated in ([4]), any 
witness approximates A to within a factor of two. This fact indicates another useful property of the Triangle 
Algorithm. 

One of the best known algorithms for determining the distance between two convex polyposes is Gilbert's 
algorithm, [H]. The connections and equivalence of Gilbert's algorithm and Frank- Wolfe algorithm, a 
gradient descent algorithm, when applied to the minimization of a convex quadratic over a simplex is formally 
studied in Gartner and Jaggi [15] . They make use of a notion called coreset, previously studied in 10,, and 
define a notion of e-approximation which is different from our notion given in (|T]). Furthermore, from 
the description of Gilbert's algorithm in [TS] it does not follow that Gilbert's algorithm and the Triangle 
Algorithm are identical. However, there are similarities in theoretical performance of the two algorithms 
and we will discuss these next. Indeed we believe that the simplicity of the Triangle Algorithm and the new 
duality theorem that inspires the algorithm, as well as the its theoretical performance makes it distinct from 
other algorithms for the convex hull decision problem. Furthermore, these features of the Triangle Algorithm 
may encourage and inspire new applications of the algorithm and further theoretical analysis, in particular 
amortized complexity of the Triangle Algorithm. In upcoming reports we shall present some such results. 

The convex hull decision problem, and its optimization form in ([5]) can equivalently be formulated as the 
minimization of a convex quadratic function over a simplex: 

n n 

min {f (x) = d{^XiVi,pf : a: £ E}, E = G E" : ^ = 1, > 0, i = 1, . . . , n}. (6) 

1=1 i=l 

A greedy algorithm for the above optimization (as well as more general smooth convex function /(a;)), as 
described in Clarkson |10| , is the following (we have changed concave maximization to convex minimization) : 

Greedy Algorithm 

• Step I. Given x' G S, let j be the index satisfying ^^q^. ^ — min{^^^-^, j = 1, . . . , n}. 

• Step II. Compute x" = argmin{/(x' + a(ej —x')) : a £ [0, 1]}, where ej is the j-th vector 
of the standard basis. Replace x' with x", go to Step I. 



Step 1 of the Triangle Algorithm and Step I of the Greedy Algorithm (also known as sparse greedy 
approximation) have in common the fact that they select an index j so that Vj will be used as a pivot 
point. Having computed such a pivot point, Step 2 of the Triangle algorithm and Step II of the Greedy 
Algorithm simply perform a line search. However, the motivation behind the selection of the index j is very 
different. The Greedy Algorithm coincides with Frank- Wolfe algorithm and Gilbert's algorithm. The Greedy 
Algorithm is algebraically motivated, while the Triangle Algorithm is geometrically motivated. The Triangle 
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Algorithm does not need to search over ah the indices to find a pivot point Vj . In its best case it finds such 
j in one iteration over the indices by performing only 0{m) arithmetic operations. In the worst-case an 
iteration will require 0{mn) operations. The Greedy Algorithm in contrast would require 0{mn) arithmetic 
operations in every iteration. This can be seen when f{x) is written as c?(j4a;,p)^. Its gradient at a point 
would in particular require the computation of Ax. 

The Greedy Algorithm generates a sequence of vectors ccj/j) G R" where X(^k+i) has at most k nonzero 
coordinates. This is advantageous when n is very large. As will be easily verifiable this property also holds 
for the Triangle Algorithm when the initial iterate po is taken to be sparse. Another property of the Greedy 
Algorithm is that if Xi^ IS the optimal solution of ([5]), then 

f{x(^))-f{x.)<^ = 0{\), (7) 

where C/ is a constant that depends on the Hessian of /, see Clarkson 10 and Zhang [39]. 

The Triangle Algorithm generates a sequence of points p'^^.^ in conv{S) that get closer and closer to p. 
This sequence corresponds to a sequence in S, where /(a^^j,-,) — d(j>'^i^y p)'^ . It p S conv{S), /(x*) — 0. 
Given any e G (0, 1), the Triangle Algorithm in K^^ < 48e~^ iterations will generate a point p^ G conv{S) 
satisfying 

d{p[K,)^P) i? = max{c?(p, = 1, . . . ,n}. (8) 

Given an index fc, by reversing the the role of k and e and solving for e in the equation = k, we get 

e = ^ 48/fc so that we may write 

dip[^^,p) < ^ -R. (9) 

Equivalently, 

/K,))<^^0(i). (10) 

In summary, when p £ conv{S) the Triangle Algorithm works similar to the Greedy Algorithm but it may 
perform faster because it only needs to find a pivot point vj as opposed to finding the minimum of partial 
derivatives in Step I of the Greedy Algorithm. When p is not in conv{S), the Triangle Algorithm can use 
any witness to give an approximation of closest point to within a factor of two, see (jlj. We may conclude 
that the Triangle Algorithm in theory is at least as effective as the Greedy Algorithm, and possibly faster 
whether approximating p G conv{S), or estimating the distance A to within a factor of two. 

In the context support vector machines (SVM) (see [3] for applications) , Har-Peled et al. [TH] use coreset 
to give an approximation algorithm, see also Zimak '3D' . We mention these because we feel that the Triangle 
Algorithm, despite some similarities with existing algorithms or their analysis, is distinct from them. It 
is quite simple and geometrically inspired by the distance duality, a simple but new and rather surprising 
property. 

1.4 Solving LP Feasibility by the Triangle Algorithm 

The LP feasibility problem is to test if the polyhedron 

= {x G R" : Ax^b, x>0} (11) 
is nonempty, where A is an to x n real matrix. The set of recession directions of Q is the set 

Res{n) ^{d: Ad = 0, d>0, d^O}. (12) 
If A — [fli, . . . , a„], it is trivial to show 

Res{n) = ^ conv{{ai, an})- (13) 
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When i?es(ri) ~ 0, it is easy to prove f2 ^ if and only if G conv(\a\, . . . , a„, —6}). In this case, given 
any e G (0, 1) we can apply the Triangle Algorithm to get an e-approximate solution: 

n 

p' = ^ aiai - an+ib G conv{{ai, ...,an, -b}) (14) 



where 
with 



d{p,0)<eR', (15) 



R' = max{d(ai, 0), . . . , d(a„, 0), d{b, 0)}. (16) 
Setting xq = {ai/an+i, ■ ■ ■ ,a„/a„+i)"^, xq >0 and from ([T5|) we have 

ei?' 

d(^a;o,&)< . (17) 

However, a lower bound on q;„-|_i is needed in order to estimate the quality of xq as an approximate 
solution of f2. A theoretical lower bound is the quantity: 

n 

Ao = mm{d{p,0) : p £ conw({ai, . . . , a„})} = min{d(Aa;, 0) : ^0:^ = 1, a; > 0}, (18) 

A computational lower bound to Ag can be derived by applying the Triangle Algorithm itself. From 
the property stated in any witness p' £ conv{{ai, . . . , an}), i.e. satisfying d(p',ai) < d{0,ai), for all 
I = 1, . . . , n, implies 

id(p',0)< Ao<d(p',0). (19) 

We prove a sensitivity theorem that given any lower bound Aq to Aq, gives the accuracy to which we need 
to solve the convex hull problem corresponding to testing if £ conv{{ai, . . . , a„, —b}). Setting bo = d(b, 0), 
our sensitivity theorem implies that for any e < A'q/2R' , we have 

d{Axo,b)<e'R', £' = 2^1 + ^) (20) 

From (PO)) it follows that if for a given eg G (0, 1) we wish to compute xq > such that d{Axo, b) < eoR', 
it suffices to have e satisfying 

In particular, it suffices to pick e < eoAQ/4i?'. Applying these, we offer a two-phase Triangle Algorithm 
for solving the feasibility problem in LP with the assumption that ^ conv{{ai, . . . , a„}): 



Two-Phase Triangle Algorithm (A = [ai, . . . , a„], b) 

• Phase 1. Call Triangle Algorithm({ai, . . . , a„}, 0) to get a witness p' £ 
conv{{ai, . . . , an})- 

• Phase 2. Starting with p' in Phase 1, call Triangle Algorithm({ai, . . . ,a„, — 6},0). 



In particular, the number of iterations in Phase 2 of the algorithm is 0(e ) = 0(eQ (i?'/Ag) ). 

Remark 4. Another lower bound to Aq can be argued to be of the order of 2~'~^^-'^\ where L is the size of 
encoding of A and b. For such bound coming from LP-type analysis, see [53]. 
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1.5 Problem (P) and Linear Programming Algorithms 

Linear programming has found numerous practical and theoretical applications in applied mathematics, 
computer science, operation research and more. In particular, the simplex method of Dantiz is not only 
a significant algorithm for solving LP but also a theoretical tool to prove many results. Ever since the 
Klee-Minty |29| example showed exponential worst-case time complexity of the simplex method, many LP 
algorithms have been invented. The trend will most likely continue. 

Problem (P) is a very special case of the LP feasibility problem. However, in fact the general LP with in- 
teger inputs can be formulated as a homogeneous case of problem (P), i.e. p — 0. The corresponding problem 
(P), may be referred as homogeneous feasibility problem (HFP), see |22j . A classical duality corresponding to 
HEP is Gordan's theorem, a special case of the separating hyperplane theorem, easily provable from Earkas 
lemma: either G conv{S), or there exists y £ M™ such that y'^Vi > for all i = 1, . . . , n, see e.g. Chvatal [7]. 
It can be justified that both Khachiyan and Karmarkar algorithms explicitly or implicitly are designed to solve 
HEP. This is because on the one hand Karmarkar's canonical formulation of LP can easily be converted into 
an HEP. On the other hand, Khachiyan's ellipsoid algorithm solves a system of strict inequalities {Ax < b) 
whose alternative system, by Gordan's theorem is an HEP {A^y = 0, b^y + s = 0, ^ y^-l-s = 1, y > 0, s > 0)). 
For this and additional results on the connections between HEP and LP feasibility, see |19| . 

By exploring the close relationship between HEP which is equivalent to finding a nontrivial nonnegative 
zero of a quadratic form, and the diagonal matrix scaling problem, Khachiyan and Kalantari |28j have 
given a very simple path-following polynomial time algorithm for LP as well as for quasi doubly stochastic 
diagonal scaling of a positive semidefinite matrix. In particular, the algorithm can test the existence of an 
e-approximate solution of HEP in a number of arithmetic operations proportional to n^'^ and lne~^. As 
approximation schemes, all known polynomial-time algorithms for LP have a complexity that is polynomial 
in the dimension of the data and in lne~^. As is well known, an exact solution for an LP with integral input 
can be computed by rounding any approximate solution having sufficient precision. Even if the complexity of 
a polynomial-time algorithms for LP would allow solving problem (P) to within e accuracy in 0(m^nlne^^) 
arithmetic operations, the Triangle Algorithm still offers an attractive alternative when the dimensions of 
the problems are large. 

Other algorithms for LP include, Megiddo's algorithm |32| with a running time that for fixed m is linear in 
n, however, has exponential complexity in m. Since Mediggo's work a number of randomized LP algorithms 
have been devised, e.g. Dyer and Erieze [13], Clarkson |5, Seidel (3^, Sharir and Welzel [32,. Kalai |20] 
gave a randomized LP simplex method with subexponential complexity bound. Matousek, Sharir, and 
Welzl |30j proved another randomized subexponential complexity algorithm for LP. See also Motawani and 
Raghavan [331. Kelner and Spielman [26] have given the first randomized polynomial-time simplex method 
that analogous to the other known polynomial-time algorithms for linear programming has a running time 
dependent polynomially on the bit-length of the input. A history of linear programming algorithms from 
computational, geometric, and complexity points of view that includes simplex, ellipsoid, and interior-point 
methods is given in Todd [55] . 

1.6 Outline of Results 

The remaining sections of the article are as follows. In Section [2l we prove several characterization theorems 
that lead to a new duality, the distance duality. We then describe the associated geometric properties and 
problems. In Section |3l using purely geometric arguments, we analyze the worst-case reduction of the gap 
in going from a new approximation to the next. Using the results in Section [3l in Section [4] we describe the 
Triangle Algorithm and analyze its worst-case complexity. In Section [S] we apply the Triangle Algorithm to 
derive a complexity for solving an LP feasibility with the assumption of that the problem, if feasible, has no 
recession direction. In Section |6l we give the analysis of the Triangle Algorithm for solving an LP feasibility 
where it is not known whether or not it has a recession direction. This corresponds to solving a general 
LP problem. The purpose and significance of the results in Section [6] lies in the simplicity of the Triangle 
Algorithm, whether or not in practice we would actually want to solve a general LP with the suggested 
approach. 
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2 Characterizations and Applications 

Throughout the section, let S = {vi, . . . , u„} be a set of points in M™, and let p G K™. 

Theorem 1. (Characterization of Feasibility) p e conv{S) if and only if given any p' G conv{S)\{p}, 
there exists Vj € S such that d{p',Vj) > d{p,Vj). 

Proof. Suppose p e conv{S). Consider the Voronoi cell of p with respect to the two point set i.e. 
V{p) — {x & M™ : d{x,p) < d{x,p')} (see Figured]). We claim there exists vj G V{p). If not, 5 is a 
subset of V{p') = {x E M™ : d{x,p) < d{x,p')}. But since V{p') is convex it contains conv{S). But since 
V{p) n V{p') = this contradicts that p E conv{S). 

Conversely, suppose that for any p' G conv{S) \ {p} there exists Vj such that d{p',Vj) > d{p,Vj). If 
p conv{S), let p' G conv{S) be the closest point to p. Since the closest point is unique, for each i = 1, . . . , n, 
in the triangle App'vi the angle Zpp'vi must be non-acute. Hence, d(p' , Vi) < d(p, Vi) for all i, a contradiction. 
Thus, p G conv(S). □ 

The following is a convenient restatement of Theoremfl] relaxing the strict inequality, d{p' , vj) > d{p, Vj). 
It has an identical proof to that theorem. 




Figure 1: Example of cases where orthogonal bisector oi pp' does and does not separate p from conv{S). 



Theorem 2. p G conv{S) if and only if given any p' G conv{S) \ {p}, there exists Vj G S such that 
d{p',Vj) > d{p,Vj). □ 

Definition 1. We call a point p' G conv{S) a witness if d{p', Vi) < d{p, Vi), for alH = 1, . . . , n. We denote 
the set of all witnesses by Wp . 

The following is a characterization of infeasibility. 

Theorem 3. (Characterization of Infeasibility) p ^ conv{S) if and only if there exists a witness p' . 

Proof. Suppose p ^ conv{S). Then by Theorem [2] there exists p' G conv{S) such that d{p',Vi) < d(j),Vi), 
for all i — But then p' is a witness. Conversely, given that p' is a witness, Theorem [T] implies 

p ^ conv{S). □ 

Thus we may conclude the following, a new duality for the convex hull decision problem. 

Theorem 4. (Distance Duality) Precisely one of the two conditions is satisfied: 
(i) For each p' G conv{S) there exists Vj G S such that d{p',Vj) > d{p,Vj). 
(a) There exists a witness. □ 

The following is a straightforward but geometrically appealing characterization of the set of witnesses as 
the intersection of open balls and conv{S). 
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Proposition 1. Let Bi = {x € M™ : d{x,Vi) < d{p,Vi)}, i = Then Wp = conv{S) n (n'^^^Bi). 

In particular, Wp is convex and lies in the relative interior of conv{S) . □ 

Figure [5] gives a case wlien Wp is empty. Figure [3] gives several scenarios when Wp is nonempty. Next, 
we state a corollary of Theorem [T] and Theorem [5] 




Figure 2: A case with no witnesses: p € conv(S). 



Corollary 1. (Intersecting Balls Property) Consider the set of open balls Bi = {x € M™ : d{x,Vi) < 
ri}. Assume they have a common boundary point p, i.e. p € D^^^dBi = {x € : d{x,Vi) = r^}. Let 
S = {vi, . . . , Vn}. Then p e conv{S) if and only if (nf^^Bj) D conv{S) — (see Figure[2]). 

Proof. Suppose p G conv{S). Pick any point p' £ conv{S) \ {p}. Then by Theorem [1] there exists Vj such 
that d{p',Vj) > d{p,Vj). But this implies p' ^ Bj, hence (n"^ji3i) fl conv{S) = 0. 

Conversely, suppose that (fl^^iBi) Ci conv{S) = 0. Thus for each p' G conv{S) there exists vj such that 
d{p' ,Vj) > d{p,Vj). Then by Theorem[5]we have p £ conv{S). □ 

Proposition 2. (The Orthogonal Bisector Property) Suppose p' is a witness, i.e. p' G conv{S) 
satisfies d{p' , vi) < d{p, Vi), for all i — 1, . . . ,n. Then, the orthogonal bisector hyperplane of the line segment 
pp' separates p from conv{S). More specifically, let c = p — p' and 7 = ^{p^p — p'^p'). If H — {x £ R™ : 
c^x = 7}, then p e 7J+ = {a; G R" : c^x> 7} and conv{S) C F_ = {2: G R" : c^x < 7}. 

Proof. It is easy to verify that p G We claim Vi G -ff-, for each i. For each i = 1, . . . ,n we have, 

d^(p', Vi) < d'^ip, Vi). Equivalently, {p' - Vi)^{p' - v,) < {p - ViY [p - Vi). Simphfying, gives 

2(p - pfv^ < [p^p - p'^p). 

Hence S C . Since i7_ is convex, any convex combination of points in S is also in . □ 

We now give a complete characterization of the witness set. 

Theorem 5. (Characterization of Witness Set) p' G Wp if and only if the orthogonal bisector hyperplane 
of the line segment pp' separates p from conv(S). 

Proof. By Proposition[2l p' G Wp implies the orthogonal bisector hyperplane of the line segment pp' separates 
p from conv{S). Conversely, suppose for some p' G conv{S) the orthogonal bisector hyperplane of the line 
segment connecting pp' separates p from conv{S). Then, in particular we have d{p',Vi) < d{p,Vi), for all 
i — 1, . . . ,n. Hence, p' G Wp. □ 

Corollary 2. Suppose p ^ conv{S) . Let 

A = d(p, conv{S)) = min{d(p, x) : a; G conv{S)}. 
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Figure 3: Examples of nonempty witness set Wp, gray areas: p ^ conv{S). 

Then any p' € Wp gives an estimate of A to within a factor of two. More precisely, 

< A < d{p,p'). 

Proof. The inequality A < d{p,p') is obvious. The first inequality follows from Theorem [5] □ 




Figure 4: The set Wp of general witness (dotted area) is a superset of Wp (gray area). 

Before we utilize the characterization theorems proved here we wish to give a variation of Theorem [T] 
The theorem shows that the notion of a witness need not be restricted to the convex hull of S. The proof is 
identical to the proof of Theorem 1 and is omitted. 

Theorem 6. p e conv(S) if and only if given any point p' ^ p, there exists Vj £ S such that d(j)',Vj) > 
d{p,Vj). □ 

We may thus give a more general distance duality as well as definition for a witness. 

Definition 2. We call a point p' e M™ a general witness if d{p' , Vi) < d{p, Vi), for alH = 1, . . . , n. We denote 
the set of all general witnesses by Wp. 

Figure |4] depicts the set Wp which of course contains Wp as a subset. The following is a corollary of 
Theorem [2] and Theorem |6] 
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Corollary 3. (General Intersecting Balls Property) Consider the set of open balls Bi = {x E R™ : 
d{x, Vi) < Vi}. Assume p € Cif^idBi. Then p £ conv{{vi, . . . , t;„}) if and only if (Hf^iBi) — 0. 

Proof. Suppose p G conv{S). Pick any point p' ^ p. Then by Theorem [B] there exists Vj such that d(jj' , Vj) > 
d{p,Vj). But this imphes p' ^ Bi, hence (nf^^Bi) = 0. 

Conversely, suppose that (D^^iBi) — 0. In particular, for each p' £ conv{S) there exists Vj such that 
d{p',Vj) > d{p,Vj). Then by Theorem [5] we have p £ conv{S). □ 

Remark 5. The general open balls property suggests that in proving the infeasibility of p we have the 
freedom of choosing a witnesses outside of the convex hulls of S. However, algorithmically it may have no 
advantage over the Triangle Algorithm to be formally described later. 



3 Reduction of Gap and Its Worst-Case Analysis 

In what follows we state a theorem that is fundamental in the analysis of the algorithm to be described 
in the next section. It relates to one iteration of the algorithm to test if p lies in conv(S) and reveals its 
worst-case performance. In the theorem below the reader may consider p' as a given point in conv{S) and 
z) as a point in S to be used as a pivot point in order to compute a new point p" in conv{S) where the new 
gap d{p" ,p) is to be a reduction of the current gap d{p' ,p). 

Theorem 7. Let p,p' ,v be distinct points in R™. Suppose d{p, v) < d{p' , v). Let p" be the point on the line 
segment p'v that is closest to p. Assume p" ^ v. Let 5 = d{p' ,p), 5' = d{p",p), and r = d{p, v) (see Figure 
[5]) . Then, 



[r, otherwise. 

Proof. Given 5 < r, consider p' as a variable x' and the corresponding p" as x" . We will consider the 
maximum value of d{x" ,p) subject to the desired constraints. We will prove 



c2 

{d{x",p): xeM™, d{x',p)=S, d{x' ,v) > r} ^ SJ 1 ~ (23) 

This optimization problem can be stated in the two-dimensional Euclidean plane. Assume p ^ p" , and 
consider the two-dimensional plane that passes through the points p,p',v. Given that d < r, p' must lie 
inside or on the boundary of the circle of radius S centered at p, but outside or on the boundary of the circle 
of radius r centered at v, see Figure [Sj circles C, C, and C" . 

Now consider the circle of radius S centered at p, C" in Figure [5] Consider the ratio S'/r as p' ranges 
over all the points on the circumference of C" while outside or on the boundary of C . It is geometrically 
obvious and easy to argue that this ratio is maximized when p' is a point of intersection of the circles C" 
and C", denoted by x* in Figure [51 We now compute the corresponding ratio. 

Considering Figure |6l and the isosceles triangle Avpx* , let h denote the length of the bisector line from 
V to the base, and let q denote the midpoint oi p and x*. Consider the right triangles Apvq and Apx*x**. 
The angle Zvpq is identical with Zpx*x**. Hence, the two triangles are similar and we may write 



(24) 



1 ^ 1^ ^ rV" " ^ ' \ 
This proves the first inequality in ([22]). 

Next, suppose 5 > r. Figure [7| considers several possible scenarios for this case. If in the triangle Apvp' 
the angle Apvp' is acute, the line segment p'v must necessarily intersect C . This implies p" is the bisector 
of a chord in C , hence inside of C". If the angle Zpvp' is at least 7r/2, then p" will coincide with v. Hence, 
proving the second inequality in ([22)) . □ 
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Figure 5: Depiction of gaps S — d{p' ,p), 5' ~ d{p",p), when 5 < r = d{p, v). 



4 The Triangle Algorithm and Its Analysis 

In this section we describe a simple algorithm for solving problem (P). For convenience we shall refer to this 
algorithm as the Triangle Algorithm. The justification in the name lies in the fact that in each iteration the 
algorithm searches for a triangle App'vj where vj & S, p' & conv{S) \ {p}, such that d{p',Vj) > d{p,Vj). 
Given that such triangle exists, it uses vj as a pivot point to "pull" the current iterate p' closer to p to get 
a new iterate p" € conv{S). If no such a triangle exists, then by Theorem[3l p' is a witness certifying that p 
is not in conv{S). The Triangle Algorithm consists of iterating the following two steps: 



Triangle Algorithm {S — {vi, . . . , u„}, p) 

• Step 1. Given p' £ conv{S) \ {p}, check if there exists a pivot point Vj £ S, i.e. 
d{p, Vj). If no such Vj exists, then p' is a witness, stop. 


d{p',v,) > 


• Step 2. Otherwise, compute the step-size 




ip-p')Tiv,-p') 

(y, — . 

d'^{vj,p') 


(25) 


Let the iterate be defined as 




„ j {I - a)p' + avj , ifae[0,l]; 
1 Vj , otherwise. 


(26) 


Replace p' with p", go to Step 1. 





By an easy calculation that shift p' to the origin, it follows that the point p" in Step 2 is the closest 
point to p on the line p'vj. Since p" is a convex combination of p' and Vj it will remain in conv{S). The 
algorithm replaces p' with p" and repeats the above iterative step. Note that a pivot point Vj may or may 
not be a vertex of conv(S). In the following we state some basic properties of the algorithm to be used in 
the analysis of its complexity. 

Proposition 3. The algorithm satisfies the following properties: 

(i) In each iteration of Step 2 if p" ^ Vj, the step-size a lies in (0, 1). 
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Figure 6: The worst-case senario for the gap S' = S*, when S < r. 



(a) Given an explicit representation of p' as a convex combination ofvi's, 

n n 

p'^^QfiUj, ^Qfj^l, Q!i > 0, Vi, (27) 

1=1 i=l 

p" can also be explicitly written as a convex combination of Vi 's, 

n 

p" - = (1 " + = " ^ (2^) 

4=1 

(^iiij iSac/i iteration of the algorithm take at most 0{mn) arithmetic operations and comparisons. 

(iv) Each Vi can be selected as an iterate p" at most once. 

(v) If Vj is selected as a pivot point more than once, then except possibly for its first selection as an 
iterate, in any subsequent selection of Vj as a pivot point the iterate pk will satisfy d{p,Pk) < d{p, Vj). 

Proof. To prove (i), note that we may have p" = Vj, occurring when the Une pv and p'v are orthogonal or 
when they make an obtuse angle (see Figure [71 rightmost case of p'). Since p' ^ p, p" cannot coincide with p' . 
Otherwise, this would imply in the triangle App'v the angle Zpp'v is non-acute, implying d{p, Vj) > d{p' , Vj). 
But this contradicts that Vj is a pivot point. Given that p" ^ Vj, in the triangle Ap'pvj the angle /-p'pvj is 
obtuse. Therefor, if p" 7^ Wj, < a < 1. 

The proof of (ii) is straightforward from the fact that < a < 1 since this implies /3i > for all i and 
they sum up to one. 

The proof of (iii) is obvious since in the worst-case in each iteration we need to compute and compare 
d(p, Vi) and d{p' , Vi) for i = 1, . . . n. 

The proof of (iv) follows from the fact that the sequence of distances d{p,pk) is a decreasing sequence. 

To prove (v), assume Vj is selected as a pivot and d{p,pk) > d{p,Vj). If Vj becomes an iterate, then the 
proof follows from (iv). Otherwise, from Theorem[71 if (i(pfc,Uj) > d{p,Vj), then d{p,pk+i) < d{p,Vj). Then, 
the monotonicity property of the sequence of distances implies the desired result. □ 

We are bow ready to analyze the complexity of the algorithm. Set 

i? = max {(i(p, Ui), i = 1, . . . , n}. (29) 
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Corollary 4. Let p,p' ,p" ,v be as in Theorem^ and r = d{p,v), S — d(p,p') , 5' — d(p,p"). If p" ^ v and 

S < r, then 

^'<^/^<^exp(-i;,). (30) 

Proof. The first inequality follows from Theorem [7] and the definition of R. To prove the next inequality, we 
use that 1 + a; < exp(a;), and set x ~ □ 




Figure 7: Depiction of gaps 5 = d{p' ,p), 5' = d{p" ,p), when 5 > r ~ d{p, v). 



Remark 6. The analysis of complexity of the Triangle Algorithm will make repeated use of ([50)1 in Corollary 
m By part (v) of Proposition [3l the number of iterations where an element Vj G 5* is selected as an iterate 
is at most n. Therefore, except for n iterations, in any other iteration of the algorithm the inequality S < r 
is satisfied so that pop will apply. Hence, for the sake of simplicity in the forgoing complexity analysis we 
will exclude the occurrence of these exceptional iterations. 

Lemma 1. Assume p is in conv{S). Pick po £ conv{S) \ {p}, let 5o — d{pQ,p). Let k = k{So) be the 
maximum number of iterations of the Triangle Algorithm in order to compute a point pk in conv{S) so that 
if Sj — d{pj , p) for j = 1, . . . ,k, we have 

Sk<^<S,, (31) 

Then, k satisfies 

k^k{5o)<\No^, No = N{5o) = {32\n2)^ <23^ (32) 

Proof. For each j — 1, . . . , k — 1, the repeated application of ([30|) in Corollary 21 the fact that for each such 
j, 6j > 5q/2, and the monotonicity of the exponential function implies 

< ^.-1 exp < exp (-^) < exp {^^)<- (33) 
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It follows that 



4-1 < (5oexp - 



32i?2 



Thus from jHO]) and ^ we have 

5k < 4-1 exp 
To have 4 < 4/2, it suffices to satisfy 



8i?2 



< 4 exp 



k6, 



32i?2 J ■ 



exp 



2 X 1 

< 



32i?2 y - 2' 

Solving for k in the above inequality implies the claimed bound in p2p . 



(34) 

(35) 

(36) 
□ 



Theorem 8. T/ie Triangle Algorithm correctly solves problem (P) as follows: 

(i) Suppose p € conv{S). Given e > 0, the number of iterations Kg to compute a point p^ in conv{S) so 
that d{p,pc) < ed{p,Vi), for some Vi £ S satisfies 



(ii) Suppose p ^ conv{S). If A denotes the distance from p to conv{S), i.e. 

A — min|(i(x,p) : x G conv{S)^ ^ 



(37) 



(38) 



the number of iterations Ka to compute a witness, a point pA in conv{S) so that d{pA, Vi) < d{p, Vi) for all 
Vi £ S, satisfies 

^8i?2 _ /24' " 



hW ^ 



(39) 



Proof. From Lemma [Hand definition of fe(4) (see (|32|)). in order to half the initial gap from 4 to 4/2, 
in the worst-case the Triangle Algorithm requires k{do) iterations. Then, in order to reduce the gap from 
4/2 to 4/4 it requires at most fc(4/2) iterations, and so on. From l\'62\i . for each nonnegative integer r the 
worst-case number of iterations to reduce a gap from 4/2'^ to 4/2'"^^ is given by 



Therefore, if t is the smallest index such that 4/2* < eR, i.e. 



< 



5o_ 
Re 



then the total number of iterations of the algorithm, K^, to test if condition (i) is valid satisfies: 
< [iVol (1 + 22 + 2^ + . . . + 22(*-i)) < [iVol -^r- ^ x 2^^'-'^ < {No + 1)-^ 



From 



we get 



Ke < ( 23:^ + 1 



251 



= 23 



5l\ 2 

i?2 



(40) 



(41) 



(42) 



(43) 



Since p £ conv(S) and from the definition of R (see 
claimed bound on in (l37l). 



we have 5o = d{p,po) < R, hence we get the 
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Suppose p ^ conv{S). With po G conv{S), assume that the algorithm generates pj G conv[S) \ Wp. Since 
p ^ conv{S), for j = 0, . . . , fc we have 

A<6,^dip,,p). (44) 
From the repeated apphcation of (15(1)) in CoroUary 2] we have 

S, < 5.-1 exp (^-^ j < 4-1 exp (^-^ j < 4-2 exp (^-^ j < • • • (45) 

Thus 

4 < Jo exp (^-—j. (46) 

To determine i^A it suffices to solve 

...p(-iJ).t. m 

Solving for k = gives the claimed bound in ([5^ □ 



5 Solving LP Feasibility Via The Triangle Algorithm 

Consider the LP feasibility problem, testing if D, is nonempty, where 

n = {xeM." : Ax = b, x> 0}, (48) 

with A — [ai, a2, . . . , a„], a.; € M™. 

It is well known that through linear programming duality, the general LP problem can be reduced to 
a single LP feasibility problem. In this section we show how the Triangle Algorithm for problem (P) can 
be modified to either prove that Q is empty, or to compute an approximate feasible point when the set of 
recession directions of Q is empty, where 

Res{n) = {d: Ad = 0, d>0, d^O}. (49) 

Definition 3. Given e, we shall say xq G M" is a e- approximate solution (or feasible point) of ft if the point 
Axq lies in Cone({ai, . . . ,a„}), i.e. 

Axq e {Ax : X >0}, 

satisfying 

d{Axo,b) <eR', (50) 

where 

R' = max{d(ai, 0), . . . , d(a„, 0), d{b, 0)}. (51) 

The following is easy to show. 

Proposition 4. Suppose ^ conv{{ai, . . . , a„}) — {Ax : Xi = l,x > 0}. Then fl = ^ if and only if 

G coni;({ai, . . . , a„, — 6}). □ 

Remark 7. It is easy to see that i?es(f2) = if and only if ^ conu({ai, . . . , a„}). In particular, if 
Res{fl) = 0, ri is a bounded set, possibly empty. Thus, in this case LP feasibility reduces to a single 
convex hull decision problem. The following sensitivity theorem establish the needed accuracy to which an 
approximate solution in conv{{ai, . . . , a„, —b}) should be computed. 
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Theorem 9. (Sensitivity Theorem) Suppose ^ conv{{ai, . . . ,a„}). Let 

n 

Ao = min{rf(p, 0) : p G conv{{ai, . . . , On})} = va\n{d{Ax, 0) : a;^ = 1, a; > 0}, (52) 



=1 



Let Ag 6e anj/ number such that < Ag < Aq. Let bo — d{b, 0) and R' as in \51\l. Suppose e > satisfies 

Suppose we have computed 

p' = (aioi H 1- a„a„) - a„+i5 G cont)({ai, . . . ,a„, -5}) (54) 

satisfying 

d{p',0)<eR' (55) 

Lei 

xo = (^,...,^)^ (56) 



Then, xq > 0, and if 



ai 


an 


[ ' ■ 

Ctn+l 


' an+1 


= 2(^1^ 


bo \ 



(57) 



d{Axo,b) <e'R', (58) 
i.e. a;o is an e' -approximate feasible point ofVL. 
Proof. Letting q' =p'/a„+i, from (IMl) and ([55)) we have 

q' = ^— = Axo ~ b. (59) 

From ([55)) we have, 

0) = d{Axo,b) < (60) 

We wish to compute a lower bound on an+i- Let 

q=- ai H h- a„. (61) 

1 - a„+i 1 - a„+i 

Note that q £ conv{{ai, . . . , a„}). Then, by definition of Aq we have, 

d(g,0) > Ao > A[,. (62) 

Let 

q = = g - 6. (63) 

1 — an+i 1 - a„+i 

From ([55]) we also have 

d{q",0)< ^ . (64) 

Applying the triangle inequality, d{u, 0) — (i(w, 0) < w), to and then using the bound in we have 

d{q, 0) - ^ 0) < d(g", 0) < . (65) 
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From dnSI) and and that d{b, 0) = bo we get 

a, 

1 - ttn + l " ' 1 - a„+i 



Aq < 6o + -j : ■ (66) 



Equivalently, 

Ao-ei?' <a«+i(A^ + M- (67) 
From ([67l) and the assumption in ([SS)) we get 

A' - ei?' a;, 

Substituting the lower bound in for q;„+i into ([SO)) imphes the claimed error bound in (1551) . □ 

Remark 8. From the above theorem it follows that the LP feasibility problem when has no recession 
direction reduces to problem (P). To get an ep-approximate solution of fl with a prescribed accuracy eo < 
Aq/2R', from (157)) it suffices to apply the Triangle Algorithm using as tolerance eoAg/2(AQ + bo). In theory, 
to estimate the number of needed iterations we need to have an estimate Aq, a lower bound to Aq. However, 
in practice given a prescribed accuracy eo we merely need to run the Triangle Algorithm until either we have 
computed an eo-approximate solution of $7, or a witness proving that it is empty. Despite this alternative, 
there is a way to get an estimate of Aq as described in the next remark. 

Remark 9. Since Aq is unknown, in practice we can use an estimate. One possible approach is first to try 
to test if lies in conv{{ai, . . . ,a„}). Assuming that fl has no recession direction, is not in this convex 
hull. Thus by applying the Triangle Algorithm we will get a witness p' such that d{p',ai) < d{p',0) for all 
i = 1, . . . ,n. Such a point p' by Corollary [2] will necessary satisfy: 

^d{p',0)<Ao<d{p',0). 

Thus by solving this auxiliary convex hull decision problem we get a witness p' whose norm can be used as 
Aq. Next we use p' as the starting iterate as it already lies in the conv{{ai, . . . , a„, —b}). 

Following the above we offer a two-phase Triangle Algorithm for solving the feasibility problem in LP 
with the assumption that ^ conv{{ai, . . . , a„}): 



Two-Phase Triangle Algorithm {A ^ [ai, . . . , a„], b) 

• Phase 1. Call Triangle Algorithm({ai, . . . , a„}, 0) to get a witness p' G 
conv{{ai, . . . , an})- 

• Phase 2. Starting with p' in Phase 1, call Triangle Algorithm({ai, . . . ,a„, —b},0). 



Theorem 10. Suppose Q is nonempty and Res(fl) is empty. Then, given any eq G (0, 1), in order to 
compute an e^- approximate solution ofQ (i.e. a solution xq > such that d{AxQ, b) < e^R' ) it suffices to set 
Aq = Q.bd{p' where p' is the witness computed in Phase 1 of the Two-Phase Triangle Algorithm. Then 
in Phase 2 of the algorithm it suffices to compute a point p' G conv{{a\, . . . , a„, —6} so that 

d{p',0) < ei?', 

where e satisfies 

^ A^, . r 1 e„ 
e < mm < — , 



2 i i?' ' (A[, + bo) 
In particular, since Aq < R' , it .suffices to pick 

A' 



4R''°- 

Then the number of iterations in Phase 2 of the algorithm, K^, satisfies < 48e^^ ~ 0(e^^) = 0(eQ ^(i^'/Ap)^). 

□ 
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6 Solving The General LP Feasibility 



Here again we consider the problem of testing if 17 = {x G R" : Ax = b, x > 0} is nonempty and if 
nonempty, computing an e-approximate solution (see ([3])). Whether or not G conv{{ai, . . . , a„}), it is well 
known that to test the feasibility of one can safely add the constraint — where M is a large 

enough constant. For integer inputs, such M can be computed, dependent on the size of encoding of A, b, 
see e.g. Schrijver [3S]. In fact it can be shown that M can be taken to be 0(2*^'^^)), where L is the size of 
ri, dependent on m, n and logarithm of the absolute value of the largest entry of A or b, see e.g. [23] . 

Having such a number M at hand, by adding a slack variable, x„+i to this constraint, testing the 
feasibility of ft is equivalent to testing the feasibility of 

n+l 

r^M = {(a;,x„+i) G R"+^ : Ax = b, ^ x, = M, x > 0, x„+i > O}. (69) 

1=1 

Dividing the equations in VIm by M and setting yi = Xi/AI , j = l,...,n + l, we conclude that testing the 
feasibility of r2j\/ is equivalent to testing the feasibility of 

h "^^ 

i=l 

Next, testing the feasibility of f^M is a convex hull decision problem we call the augmented (P): testing if 
p G conv{S), where 

P=j^, 5"= {ai,...,a«,a„+i}, a„+i = 0. (71) 

We may solve the augmented (P) in two different ways, directly by solving a single convex hull decision 
problem, or by solving a sequence of such problems. We will describe the two approaches and analyze their 
complexities. First, we prove an auxiliary lemma, an intuitively simple geometric result. 

Lemma 2. Given u^w € M™, we have: 

sup|c?(u, — ): /i G [1, oo)} = max |o?(w, t/j), (i(w, 0)}. (72) 

Proof. Consider the maximum of the function f{t) — (P{u,tw) over the interval t G [0,1]. Since f{t) is 
convex, its maximum is attained at an endpoint of the interval. Letting t = l//i, the proof is complete. □ 

Theorem 11. Suppose that a positive number M satisfies the property that f2 is feasible if and only if ^Im 
is feasible. Then, solving the augmented (P) to within accuracy of ejM is equivalent to testing the feasibility 
of fl to within accuracy of e. More specifically, by applying the Triangle Algorithm to solve the augmented 
(P) (see (|7ip). in OimnM^e^^) arithmetic operations, either we compute a witness p' G conv{S) proving 
that b/M ^ S (hence — $), or a point p' G conv{S) such that for some i = l,...,n + l we have, 

d{b,Mp') < €m&x{d{b,a,),d{a^,0)} < 2R'e, (73) 

where R' is as in i51\) . Equivalently, p' = Ayo, for some yo > 0, and if xq — Myo, then 

d{b,Axo) < emax{d(6,a,),d(a,,0)} < 2R'e. (74) 

Proof. If solving the augmented (P) to accuracy e/M leads to a witness p' G conv{S) proving that b/M ^ S, 
then ri is empty. Otherwise, by Theorem [8] the algorithm leads to a point p' G conv{S) such that for some 
? = l,...,n + lwe have 

d(A p') < _Ld( — ,a,). (75) 

Multiplying the above by M, applying Lemma [5J and since Md{b/M,p') — d{b,Mp'), we get the proof 
of the first claimed inequalities in ([75)) and (|74p . The bound 2i?'e follows from the triangle inequality, 
d{b,a,) < d{b,0) + d{ai,0). □ 
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Instead of solving the augmented (P) with an a priori estimate for M, we may solve it as a sequence 
of augmented (P)'s with increasing estimates of M that successively doubles in value. To describe this 
approach, first consider the following. 

Given > 0, let — Select any p' G conv{S) as the iterate and let 

fia=snp{^: d{pf_,, at) > d{p' , ai), Vi = 1, . . . , n + l}. (76) 

Clearly, < < oo- According to Theorem [3J and the definition of /xg, for any /i £ (0, /io), p' is a witness 
proving that p^ ^ conv{S). From Lemma [51 we know /iq ^ 1- Given e > 0, Set /i = /xq and consider the 
following iterative step: 



Iterative Step. Let — e/fi. Use the Triangle Algorithm to solve the augmented (P) with 
p^ — b/ fi to either compute p'^ G conv{S) such that 

d{p^i,p'^) < £p.'aiax{d{p^,at),d{at,0)}, for some i, 

or a witness p'^ G conv{S) proving that ^ S. In the latter case replace fi with 2/i and repeat. 



Theorem 12. Repeating the iterative step, either we compute an e- approximate feasible point of fl, or 
a witness to the infeasibility of the augmented (P) with fj, > M . In the latter case is empty. More 
specifically, if the algorithm requires r iterations of the iterative step, its arithmetic complexity is bounded 
above by (641n2)mne~^2^''i?'^ = 0{mn2'^^e^^) . Furthermore, r < [logj M] . 

Proof. Since fio > 1, the number of augmented (P)'s to be solved is bounded by i = [log2M]. With the 
initial value of fi — we solve the augmented (P) to either compute an approximate solution in $7 to within 
accuracy of e, or a witness to the infeasibility of il^. By Theorem |5] this takes 0(e~^) arithmetic operations. 
If f2^ is infeasible, we double /i and test the feasibility of new il^. Then by (j32p in Lemma[I] the number of 
iterations needed to halve the current error that is known to exceed e is bounded by 

Tf - 

(32 In 2)^, i? = niax{max{d(6,a,),d(aj,0)} : i = 1, . . . ,n + 1} < 2R' , (77) 

where the bound 2R' is from the fact that d{b, Oi) < d{b, 0) + d{ai, 0), Repeating this r times we get the total 
complexity bounded by 

2 2 /2 

K, < (321n2)4-(l + 2^ + 24 + --- + 22n < (32 In 2)5^22''+i = (64 In 2)^22'^+! =0 

□ 

Concluding Remarks and Future Work. In this article we have described several novel characteri- 
zation theorems and a very simple algorithm for the convex hull decision problem, the Triangle Algorithm. 
The Triangle Algorithm is quite straightforward to implement and its main step is the computation of dis- 
tances and their comparisons in identifying a pivot point or a witness. The simplicity of the algorithm, 
its theoretical properties, and the distance duality will lead to new applications. Already in |21j we give a 
version of the Triangle Algorithm for the case of problem (P) where we wish to locate the coordinates of a 
point p having prescribed distances to the sites S = {vi, . . . , We call this the ambiguous convex hull 
problem since not only the coordinates of the point are unknown, so is the existence of such point. Despite 
this ambiguity, a variation of the Triangle Algorithm can solve the problem in 0{mne^^ Ine^^) arithmetic 
operations. We now make several comments. 

(i) While the Triangle Algorithm works with distances, it requires only the four elementary operations 
and comparisons, no square-root operation is required. Note that the constant in the worst-case analysis in 
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iteration complexity in (I37p is an absolute constant and is independent of the initial point pq. This is due 
to the fact that the Triangle Algorithm seeks to reduce the relative error, d{pk,p)/d{p,Vj), a more sensible 
measure for the problem than seeking to reduce the absolute error, d{pk,p)- 

(ii) The iteration complexity estimate of 0{e~^) is under the pessimistic assumption that in every iteration 
it has the worst-case performance, see (P^ . the relationship between S and S' in Theorem[71 Hence, in practice 
a much more efhcient complexity should be expected. 

(iii) In the worst-case, each iteration requires 0{mn) arithmetic operations. However, in the best case, 
when p € conv{S), the algorithm can find a pivot point by searching a constant number of ViS, taking 
0{m) arithmetic operations. There are some obvious options and variations of the algorithm. For instance, 
given the current iterate pk^ we can use the pivot point Vj to correspond to the smallest index j satisfying 
d{pk, Vj) > d{p, Vj). If no such j exists, pk is a witness to the infeasibility of p. Alternatively, we can choose 
among all candidate pivot points Vj, the one that minimizes the gap, the distance between the corresponding 
p" and p. This optimizes on reducing the gap but at the cost of more computation. This type of approach is 
similar to a more general version of Frank- Wolfe algorithm where all the indices are examined, see Clarkson 
[10) . Algorithm 1.2. However, in our case we suggest using only the indices that correspond to a pivot point 
(i.e. d{p',Vj) > d{p,Vj)). There are in fact many options for selecting a pivot, including relaxations in the 
computation of the closest point along the line segment p'vj , or stronger versions where we choose a direction 
in the span of directions from p' to all candidate pivot points. There are many worthy experimental and 
theoretical considerations and will be considered in a separate article. The initial iterate po can be taken 
as any of the Vi's, or the center of conv{S), namely Yl7=i '^il''^- The latter choice would of course not be a 
sparse approximation. 

(iv) Analogous to polynomial-time algorithms for linear programming with integer inputs, given that 
p S conv{S), after d{pk,p) has been sufficiently reduced, pk can be rounded into an exact solution. This 
means the a^'s in the representation of pk as a convex combination of Vi& can be rounded to an exact value 
in representing p as a convex combination of ViS. 

(v) Given an iterate pk , rather than using any representation of it as a convex combination of points in S 
as described in (I28p . we can use a representation as a convex combination of at most m-l- 1 of the ViS. From 
LP theory this is possible since each point in conv{S) can be represented as such. Such representation for 
Pk, analogous to a basis representation in the simplex method, when d[pk,p) is sufhciently small may allow 
identifying a subset of m -I- 1 points where p can be written as a convex combination of. 

(vi) It may be possible to define a path- following version of the Triangle Algorithm for solving (P) as 
follows: to test if p G conv{S), we consider p -\- tu, where u ^ is a given auxiliary point in R™, and t is 
a positive real parameter, initially selected to be large enough so that p + tu is infeasible to conv{S) via a 
witness. The parameter t is then successively reduced to zero. 

(vii) The Triangle Algorithm for (P), when applied to each Vi with Si = S — {vi} gives an approximation 
algorithm for solving the irredundancy problem. 

(viii) The Triangle Theorem can also be considered as an algorithm that solves the intersecting open 
balls problem, testing if a given set of open balls in R™ that are known to have a common point on their 
boundary, have a nonempty intersection. 

Finally, the results in this article lead into new research problems and research projects, including compu- 
tational testing, average case analysis, amortized complexity, generalization of the characterization theorems 
and the Triangle Algorithm, as well as its specialization to linear programming problems with special struc- 
tures, and combinatorial and graph optimization problems. As generalizations, it is possible to state (P) and 
the corresponding characterizations when the set S consists of one or a finite number of compact subsets of 
M™. For instance, polytopes, balls, or more general convex bodies. Moreover, (P) can also be defined over 
other cones. For instance, a canonical problem in semidefinite programming described in [24j is an analogue 
of (P) over the cone of positive semidefinite matrices. We hope to investigate these in future work. 
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